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Abstract 

The present paper numerically analyzes a passive cooling system using enclosures with different 
geometries filled with thermal conductivity-enhanced phase change material (PCM). A numer- 
ical code is developed using an unstructured finite-volume method and an enthalpy-porosity 
technique to solve for natural convection coupled to a solid-liquid phase change. Five geome- 
tries containing the same volume of PCM are compared while cooling the same surface. The 
unsteady evolution of the melting front and the velocity and temperature fields is detailed. Other 
indicators of cooling efficiency are monitored, including the maximum temperature reached at 
the cooled surface. The computational results show the high impact of varying geometry: a 
maximum temperature difference as high as 40 °C is observed between two of the enclosures. 
The best efficiency is obtained for an enclosure shifted vertically relative to the cooled surface. 
Other findings and recommendations are made for the design of PCM-filled enclosures. 

Keywords: Passive cooling, phase change material, melting, natural convection, 
enthalpy-porosity method, unstructured finite-volume method. 

1 Introduction 

The melting of phase change materials (PCM) coupled to natural convection in enclosures 
has been studied extensively. This situation is encountered in many technical applications, 
such as latent heat storage systems [TJ [2] and thermal insulation for buildings [31 H]- 
The melting of PCM is also used to control the temperature of the surface of electronic 
components that release instantaneous or periodic high density heat fluxes to moderate 
the need for classical cooling devices. There has been increasing interest in this type of 
passive cooling for electronic circuits, such as chipsets, laptop processors or graphics cards 
[5, 6J. These elements are continuously becoming smaller, and their released heat densities 
are increasing; therefore, they require efficient, economic, and silent cooling systems. A 
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particular application for such systems is the cooling of electrical devices inside shelters 
located in difficult-to-access sites and subject to periodic temperature changes. In these 
situations, cooling solutions that do not use moving mechanical elements, such as fans, 
are more suitable to avoid the need for frequent maintenance. 

Several authors have studied phase changes and natural-convection-dominated melting 
in enclosures. Beyond the development and application of numerical models and experi- 
mental techniques, these authors have proposed several methods and concepts to enhance 
the heat transfer and melting rate. Gong et al. [7j demonstrated the positive effect of 
inverting a square PCM container heated from the side during the melting process when 
thermal stratification in the melt slows the melting rate. For the same geometry, Kho- 
dadadi and Hosseinizadeh [8] numerically investigated the effect of adding nanoparticles 
to the PCM. They showed that when the enclosure is cooled laterally, the rate of freezing 
increases because of the presence of nanoparticles. Starting from a study of PCM melting 
in a square enclosure heated from below, Fteiti and Ben Nasrallah [5] examined the im- 
pact of the aspect ratio of the PCM- filled enclosure and found that flat enclosures exhibit 
faster melting but have a lower asymptotic limit of the total melt ratio (they imposed a 
temperature boundary condition equal to the melting temperature at the top side of the 
enclosure and a lower temperature at the bottom). Hernandez- Guerrero et al. [10] also 
studied the impact of the aspect ratio with a different numerical model and found similar 
qualitative results for the case of tall enclosures. It is worth noting that the two latter 
works compared differently shaped enclosures that contained varying quantities of PCM. 

Several researchers have considered the cooling of a vertical surface that dissipates heat 
flux using a tall enclosure filled with PCM. Binet and Lacroix [11] performed numerical 
and experimental studies to analyze the impact of the positions of three heat sources, 
their sizes and the aspect ratio of the enclosure itself. A similar numerical study was 
conducted by Krishnan and Garimella [12] on the conjugate heat transfer through the 
enclosure walls. Pal and Joshi [13] studied a uniformly dissipating heat source using 
experiments and numerical modeling, and they established correlations for the temporal 
evolution of the parietal heat transfer and the melt fraction. Huang et al. [H] used a 
three dimensional (3D) numerical model to investigate the cooling of photovoltaic cells 
by a PCM enclosure equipped with metallic internal fins. 

A large amount of research has been conducted to investigate the applicability of 
cooling by PCM for electronic and electric devices. For example, the work in [15] focused 
on the cooling of a mobile device using embedded PCM enclosures, [TBJ [TTJ HBJ E] focused 
on PCM-based heat sinks, while in [JJJ] PCM was used to control the temperature of 
Li-ion batteries. 

Unlike for thermal insulation, the efficiency of passive cooling with PCM is closely 
related to the rapidity of its melting: the temperature of the heat-dissipating surface does 
not rise rapidly because of the absorption of the latent heat of fusion. 

The speed with which the PCM melts in a laterally heated enclosure is controlled by 
the conductivity of the solid PCM during the early stages and by the natural convection 
currents in the melt during the later stages. This natural convection flow depends on 
the thermo-physical properties of the PCM and the geometry of the enclosure. Thus, the 
shape and extent of the enclosure containing the PCM have an important effect on the 
kinetics of the melting front and, therefore, on the temperature of the cooled wall. 

In this work, we are interested in the cooling of a vertical surface releasing a fixed 
heat flux by means of a PCM-filled enclosure that entirely covers the surface. The aim 
of this study is to propose and examine different geometric shapes and relative positions 
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of the enclosure. Numerical simulation is a well-suited tool for this purpose, especially 
in the early stages of design. It allows different test situations to be compared before an 
experimental prototype is constructed. A numerical model based on a fixed grid method 
and applied to an unstructured finite-volume formulation has been developed for this 
purpose. 

In the following sections, the studied system is described with the different considered 
geometries and the flow conditions are then detailed. Next, the physical model and 
the numerical method are presented, and the implemented code is validated. Afterwards, 
detailed results are presented concerning the evolution of the velocity and the temperature 
fields, the melting front and the rate of fusion. Finally, the geometric impact on the 
coupled effects of natural convection and phase change is analyzed. Conclusions are then 
drawn based on these results. 

2 The studied system 

In this study, we are interested in the case of a vertical wall releasing a uniform heat flux 
at a high rate but for a limited duration. We aim to maintain this surface (designated 
Sh) at a sufficiently low temperature to avoid damaging the electrical component behind 
it. To achieve this objective, we provide this surface with an enclosure filled with a PCM. 
The surface to be cooled constitutes one face, or a part of a face, of this enclosure while the 
other faces are exposed to the external air and are subject to natural convection. When 
the heat release from the surface Sh begins, its temperature rises until it reaches the 
fusion temperature of the PCM, which causes the material to melt. Natural convection 
currents take place progressively in the melt, and the generated flow ensures that heat 
is transferred from the surface Sh to the solid PCM and maintains the melting process. 
The strength of the fluid flow involved in this process is closely related to the shape of 
the formed melted area, the extent of which is limited by the boundaries of the enclosure. 
Therefore, we expect that the efficiency of the heat transfer and the rate at which PCM 
melts, strongly depend on the shape of the enclosure. This simple passive cooling device is 
appropriate for intermittent heat release and is studied here for a uniform heat flux rate of 
10 4 W.m~ 2 . For continuous heat release and higher heat flux densities, more sophisticated 
solutions could be considered such as heat pipes [20J. 

The "natural" or "intuitive" choice of the shape of the PCM container can be re- 
garded as the rectangular geometry which has been extensively studied in the literature 
for its simplicity. It can also be considered that the circular shape has the advantage of 
the absence of corners that may retain non melted PCM. Furthermore, the "intuitive" 
choice of the relative position of the container to the heating surface would be a centered 
position. Thus, in this study we aim to verify the relevance of these different choices. 
Consequently, we proceeded to investigate the impact of different enclosure shapes on 
the natural convection flow coupled to the solid-liquid phase change. The five different 
geometries studied contained the same volume of PCM. Fig. [1] shows the computational 
meshes used to model the enclosures (the chosen meshes densities are described in section 
I4.4p . In this figure, the location of the surface Sh is indicated in gray. The area of Sh is 
the same for all of the enclosures. These geometries (detailed in the next section) were 
not all considered at the beginning of this work but were designed progressively given the 
results obtained from the first tested shapes. However, the results obtained for the five 
enclosures will be presented together. Some of the enclosures are rectangular, whereas 
others have rounded corners. Barlett et al. [2T] have studied natural convection in both 
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rectangular and rounded-corner enclosures and have shown that the latter shape substan- 
tially improves heat transfer. Thus, we considered rounded shapes in the present study 
and investigated their potential role in situations involving a phase change. 

A known drawback of the considered cooling system is the relatively low thermal con- 
ductivity of common PCMs (e.g., wax paraffin), which may cause a high temperature 
peak at the surface in the case of high heat flux early-on in the heat release, when the 
heat transfer occurs only by conduction. To overcome this limitation, several solutions 
have been investigated in the literature, such as utilizing fins [B], including a high ther- 
mal conductivity metallic foam [BJ 122] or adding nanoparticles. Extensive reviews of all 
these solutions can be found in J23J El]- In our study, we consider a solution that in- 
cludes graphite nanoparticles with high thermal conductivity in the PCM. These types 
of materials have experienced extensive development and exhibit high potential for use 
in several applications. Many researchers [2"5"| I2"6] have demonstrated that the thermal 
conductivity of these products can reach unity (SI units). As indicated above, in the 
present study, we focus on the impact of the enclosure shapes on the melting kinematics, 
natural convection strength and temperature extremes at the cooled surface. Therefore, 
we choose not to consider the details of the mixture of PCM and nanoparticles and the 
relationship between the percentage of nanoparticles and the thermophysical properties 
of the obtained PCM, such as viscosity, effective conductivity and specific heat capacity. 
Instead we consider a model PCM with constant thermophysical properties but enhanced 
thermal conductivity (see next section for details). 

3 Flow conditions and PCM properties 

As stated above, in this work, we study the melting of PCM inside five enclosures of 
different shapes that surround the same volume. We limit our investigation to two- 
dimensional configurations that can be modeled in a vertical plane. One of the sides of 
the first two enclosures on the left in Fig [1] (a and b) exactly matches the surface Sh- 
The first enclosure is rounded (a half disc), and the second is rectangular. The three 
other enclosures (c, d and e) have sides that vertically exceeded the surface Sh, and thus, 
their widths are smaller than those of enclosures (a) and (b). Enclosure (c) is oblong 
with rounded corners, whereas enclosure (d) is rectangular. The last enclosure (e) has 
the same geometry as enclosure (d), but it is translated upward such that its left and 
right sides exceed the surface Sh only at the top. The surface Sh has a fixed height of 
H = 5 cm, whereas the dimensions of the five enclosures are adjusted to contain the same 
volume. These dimensions are given in Tab. [T]in terms of maximum width and height. 

Tab. 1: Maximum widths w max and heights h max of the five enclosures (in cm). 



Enclosure 






(a) 


2.5000 


5.0000 


(b) 


1.9635 


5.0000 


(c) 


1.5742 


6.5742 


(d) 


1.4933 


6.5742 


(e) 


1.4933 


6.5742 



The rate of heat flux at the surface Sh is fixed to a constant value, q" = 10 4 W.m~ 2 , for 
all cases. The heat transfer to the external air at the other boundaries is considered to oc- 
cur by natural convection and is modeled by a heat transfer coefficient h = 30 W.K~~ 1 .m~ 2 , 
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Fig. 1: The geometries of the five studied enclosures: a,b, c, d, and e (from left to right) and 
corresponding meshes (in gray: the position of the surface S^)- 



whereas the external air temperature is assumed to be equal to the initial temperature of 
the PCM, = T = 20 °C. 

We consider a model PCM in this study composed of a pure paraffin wax and graphite 
nanoparticles. To simplify the analysis in the study, the properties of this PCM are 
considered to be identical in the solid and liquid phases and non-temperature dependent 
(except in the body forces term). The values of the adopted properties are listed in Tab. 
El A maximum temperature difference AT characterizing the flow can be obtained by 
considering the final or maximum admissible temperature that reaches the surface Sh, 
which corresponds to the temperature at which an electronic device is damaged (i.e., 90 
°C): AT = 90 — 20 = 70 °C. The Rayleigh number based on this temperature difference 
and on the height H of the surface Sh is Ra = 2.7 x 10 7 , which corresponds to a laminar 
flow regime. The Stefan number corresponding to this limit situation is Ste = 0.87. 



Tab. 2: Thermophysical properties of the model enhanced-conductivity PCM 

Dynamic viscosity (/x) 5 x 10~ 3 Pa.s 

Density (p) 800 kg.m~ 3 

Thermal conductivity (k) 1 W.m^ 1 K^ 1 

Specific heat (C p ) 2500 J.kg' l K~ x 

Latent heat (L) 200 kJ.kg- 1 

Thermal dilatation coefficient (/3) 10 -3 K^ 1 

Prandtl's number (Pr) 12.5 

Melting temperature (T m ) 20 °C 
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4 Computational modeling 
4.1 The governing equations 

The unsteady equations governing the flow of incompressible non-isothermal fluids are 
solved over the entire computational domain. The three-dimensional equations of con- 
servation of momentum, mass and energy (in terms of temperature) are considered valid 
for both solid and liquid phases, which are distinguished by a liquid volume fraction / 
that takes values and 1, for solid and liquid phases, respectively These conservation 
equations are considered in their integral forms: 

/ p U dV + [ p UU ■ n dS = - [ Vp dV + [ r ■ ft dS + [ S v dV (1) 
Jv Js Jv J s Jv 

U-ndS = (2) 

s 

/ p C P T dV + [ pC p TU-ndS = [ k VT • n dS + [ pL^-dV (3) 

ut Jv JS Js Jv ut 

where U is the velocity vector, p the pressure and T the temperature, r, is the viscous 
stress tensor for a Newtonian fluid: 

t = n (VU + (VUf) (4) 



The integration occurs over a volume V surrounded by a surface S, which is oriented 
by an outward unit normal vector n. The source term in the momentum conservation 
equation (Eq. (jTD) contains two parts: 

Su = p(3(T-T ref )g + AU (5) 

where ft is the coefficient of volumetric thermal expansion and g the acceleration of gravity 
vector. The first part of this source term represents the buoyancy forces due the thermal 
dilatation. For sake of simplicity, the Boussinesq's approximation is used in this com- 
parative study that focuses on the impact of different geometries. The T re f temperature 
was chosen as the mean temperature of the PCM liquid phase and was recalculated at 
each time step. Therefore, T re f is more representative of the temperatures in the liquid 
phase throughout the whole process, especially when all of the PCM has melted. The 
second part of the source term (Eq. \5§ is a penalization term that ensures zero velocity 
in the computational cells where the PCM is solid, i.e., where / = 0. The penalization 
coefficient A, based on the Carman-Koseny relation for a porous medium, is written as a 
function of / as follows: 

where C = 1.6 x 10 6 , and e = 10~ 3 . This coefficient ensures a smooth transition between 
solid and liquid media. In the case of pure PCM, the phase change front is sharp, and 
the transition takes place over only one computational cell in the direction of the front 
displacement. 

In Eq. (jnj), the last term of the right hand side (RHS) introduces the effect of the 
latent heat of phase change L into the energy equation. This effect is accounted for by 
considering the variation in time of the liquid fraction /. 
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Fig. 2: A computational cell C and one of its neighbors Nb- 

This general three-dimensional (3D) model can be restricted to deal with 2D com- 
putations (e.g., in the (x,y) plan) without any change by considering a single layer of 
computational cells (in the z direction) and by neutralizing the top and bottom faces 
(with respect to z). It is this approach that was used for all the computations presented 
in this study, thus, the considered meshes are composed of one layer of hexahedra with 
or without prisms. 

4.2 The numerical method 

The conservation equations ((U [2] and |3]) are solved by implementing them in an in- 
house code called Tamaris. This code has a three-dimensional unstructured finite-volume 
framework that is applied to hybrid meshes. The variable values (U, p, T and /) are stored 
in cell centers in a collocated arrangement. The cell shapes can vary (e.g., tetrahedral, 
hexahedral, prismatic or pyramidal). 

To describe the discretization method used in the code, we can write Eqs. (JTJ and (J3J) 
in the generic convection-diffusion form with respect to a conserved variable 0: 

^- f P (j)dV+ f p<pU -ndS = [ TV0 • n dS + [ S$ dV, (7) 
Jv Js Js Jv 

where T is a diffusion coefficient and a source term. The spatial schemes used to 
approximate the diffusive and convective fluxes are both second-order accurate. The 
diffusion term is discretized by approximating the surface integrals with a sum over all 
cell faces / (Fig. |2j|: 

/ rV0-ndS = Y^T f A f {V$) r n f , (8) 
Js f 

where Af is the area of face /. For unstructured meshes, orthogonality is an exception, and 

it needs to be handled correctly. Thus, the normal gradient (V<p) / • Hf is decomposed into 
an implicit contribution that uses the values of </> at the centers of the two cells sharing 
the face / (the first term on the RHS of Eq. (J9])) and a non-orthogonality correction 
term treated explicitly by a deferred approach to preserve the second-order accuracy of 
the centered differencing. We use the over-relaxed decomposition suggested by [27] to 
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enhance the convergence properties of the discretization of the diffusive term: 

<ft) r « l = +^* + Wl.L-*) ( 9) 

a - n/ y a - nf J 



d is the vector joining the centers of the two adjacent cells (see Fig. [2]). The average 

gradient V0 is interpolated from the gradients of these neighboring cells. 

The gradients of the variables at the cell centers are computed by Gauss' theorem: 

v}=~ I <f>ndS=±J2 ( t ) f A fnf, (10) 
Js f 

where 0/ is the mean value of the variable interpolated using the values at the centers of 
two cells sharing face /: 



/ = ^ c +(i-O<K with £ = = (li) 

Once the gradient is calculated for all computational cells, the values are used to determine 
a new estimate of 6t as follows: 



(<K + Vjiv, • N3) + (0c + Vj c • Cf)] (12) 



These new values of 0/ are used to re-compute the gradients more accurately using Eq. 

(HDD [2E]. 

Convection terms are also transformed into a sum over faces / by decomposing the 
surface S: 

[ P <piJ-ndS = V(p0A) f U f ■ n f , (13) 
Js f 

where the face values 0/ require appropriate interpolation to be accurate and bounded. 
Thus, we use the non-linear high-resolution (HR) bounded scheme CUBISTA by Alves et 
al. [21] in the 7 formulation of Ng et al. [30J, where they expressed 0/ is a function of 
the upwind (UP) value of and its centered differencing (CD) value: 

0^ = 0^ + 7(0^-0- (14) 

The coefficient 7 is determined for each face based on the local shape of the flow solution 
using the normalized variable diagram (NVD) framework and observing the convection 
boundedness criterion (CBC) [31]. The first term of the RHS of Eq. (TH|) is accounted 
for implicitly, whereas the second term is treated explicitly with the deferred- correct ion 
practice. 

The pressure- velocity coupling is ensured by the SIMPLE algorithm [32], whereas the 
mass fluxes at the cell faces are evaluated using the Rhie-Chow interpolation [33J to avoid 
pressure checkerboarding. The implicit three-time-step Gear's scheme [31] of second-order 
accuracy is used to discretize the unsteady terms: 

ir p4>dV = 3W):-^w;-+wry (15) 

dt Jy At 
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►( New time step J 

T 

f Resolve for U , Eq. (1) Y^- 

T 

f Resolve for pressure correction J 

T 



c 



jf Correct p and IT J 

T 



Resolve for T. 



Eq. (3) 



~~T 

^ Update liquid fraction/', Eqs. (17,18) ^ 




Fig. 3: The liquid fraction updating procedure as included in the SIMPLE algorithm. 



The superscript n stands for the current time step and At 



t n -f 



n-l 



is the time step. 



The RHS of Eq. ((?]) is taken at time t n . By applying the former discretizations, the 
generic conservation Eq. ([7]) transforms into the algebraic form: 



N b 



b c . 



(16) 



Within each iteration of the SIMPLE algorithm, after the resolution of the momentum 
equation and of the Poisson equation for the pressure correction [281 E2], the energy 
equation is solved, and the fluid fraction / is updated. These last two steps are repeated 
until the variation of / is sufficiently small; then, the next SIMPLE iteration starts unless 
the convergence for U, p and T is achieved, in which case a new time step is considered. 
A diagram of this algorithm is given in Fig. |3j The resolution of the energy equation is 
integrated in the SIMPLE iteration to take into account the high level of its coupling with 
the momentum equation through the body forces term. Additionally, the liquid fraction 
correction is iterated with the energy equation resolving to tightly couple the temperature 
field with the liquid fraction field. 

In this study, the liquid fraction / is updated by the a new source' 1 algorithm proposed 
by Voller [35], where the new value of / at iteration k + 1 and in cell c is calculated as 
follows: 

At all 



f, 



k+l 



pLV 



(17) 



where T m is the melting temperature of the PCM, and is the coefficient of T c in the 
discretized Eq. ffTB"]) that corresponds to temperature. This update is followed by an 
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0.02 0.04 0.06 0.08 0.1 200 400 600 800 1000 

x [m] time [s] 

Fig. 4: Comparison of the obtained melting front positions at several instants (left) and of the 
evolution of the total liquid fraction in the enclosure (right) with the benchmark results 
of Hannoun et al. 



overshoot / undershoot correction: 



if / c fc+1 < 

1 if f* +1 > 1. 



Following the "new source" algorithm [35], the energy equation is penalized in the 
computational cells belonging to the phase change front (0 < f c < 1) to ensure that 
the temperature at these cells is equal to the melting temperature T m . This procedure is 
performed by adding a penalization source term, equal to 10 9 xT m , to the energy equations 
corresponding to these cells. This practice accelerates the convergence of the / updating 
algorithm, and the number of iterations needed to reach convergence (|/ c fc+1 — / c fc | < TOL) 
may be lowered to 1 or 2 depending on the size of the mesh and time-steps. The value 
of TOL is fixed in this work to 10~ 4 . At each iteration, the discretization technique 
presented above leads to a linear system of algebraic equations in the form of Eq. ([161) 
with a non-symmetric sparse matrix for each variable. These linear systems are solved 
using an ILU-preconditioned GMRES procedure implemented in the IML++ library |36| . 
In the scope of this work, all of the computational meshes were generated using the 
open-source software Gmsh [37] . 



4.3 The code validation 

The present code has been successfully validated in some situations involving flow and 
heat transfer, as in [38, 39J. An additional code validation concerned with pure natural 
convection is given in|AJ We focus here on validating the code when applied to the case of 
melting of a pure PCM coupled to natural convection in the melt. The chosen test case 
is the 2D numerical benchmark presented by Hannoun et al. [40J , which involves melting 
tin in a square enclosure subject on one side to a temperature higher than the melting 
temperature. The authors have presented extensive results obtained using a second-order 
accurate finite- volume method in structured meshes as fine as 600 x 600. 
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We carry out a numerical simulation in the same conditions as those described in [40] 
with the same material properties using a uniform grid 200 x 200 in size. In Fig. HJ we 
plot the melting front at different times. These lines correspond to the isolines of / = 0.5. 
The melting fronts at times 450 s and 700 s present a wave-like shape that is due to the 
presence of three and two fluid recirculations (or rolls) respectively. The top two rolls 
have merged into a single one at about t = 480 s. We also present the total liquid fraction 
Fi in the entire square enclosure, which is calculated as follows: 




where V c is the volume of a computational cell c, and the summation is over all cells in 
the computational domain. 

All these results are compared to those in [40] and exhibit satisfactory correspondence. 

4.4 Mesh size-dependence study 

To choose the mesh size with the best compromise between accuracy and computational 
cost, we conduct a mesh size- dependence study to determine the number of computational 
cells of the mesh needed to achieve satisfactory accuracy. The rectangular enclosure (b) 
undergoing an unsteady melting process in the same conditions as those described in 
section [3] is chosen as an example to explain the mesh selection process. Three meshes 
(mi, vri2 and m.3) of increasing size that contain 6000, 11660 and 24000 cells, respectively, 
are considered. The results obtained with the three meshes are compared based on the 
time evolution of the mean values and on the instantaneous spatial fields (i.e., melting 
front position, temperature and liquid fraction). 
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Fig. H^a) gives the form of the melting front obtained by the three meshes at time 
t = 200 s. These fronts are obtained by plotting the / = 0.5 isovalue line of the liquid 
fraction field. The fronts obtained by the three meshes are almost superimposed on each 
other, which indicates that the mesh size has little influence on the front position. More 
significant differences are observed in the temperature field, as seen in Fig. EJb), which 
shows the variation of the temperature along a horizontal axis passing through the center 
of the enclosure. For clarity, the plot is limited to the liquid region, where the temperature 
is not constant. However, the results of the three meshes are comparable. To achieve a 
quantitative comparison, we consider 50 points along this cutline and calculate the local 
error in the temperature at each point j as follows: e* = (Tj — Tj)/Tj x 100, where 
i = l,2 stands for meshes nil and the results are compared to those of the finest 
mesh 777.3. Thus, the mean value of the error is e 1 = Y2j e j/50, and the maximum error is 
e max = m aXj e* . They are both reported in Tab. [3J 



Tab. 3: Relative error of the local values of T along the central horizontal cutline at 
t = 200 s using the results given by the finest mesh m 3 as the baseline. 



mesh 


e 


&max 


mi 


2.2% 


8.1% 


m 2 


1.2% 


6.1% 



In addition to the spatial results, we compare the time evolution of the total liquid 
fraction in the enclosure F\ (Eq. (TT5])) and the mean temperature T sm of the surface Sh, 
which is calculated as follows: 



i^ A ' T ') 



Tsm = ( > , A f Tf ) , (20) 

where Af is the area of a mesh face / located on the surface Sh- Relative errors of meshes 
mi and m 2 compared to are calculated at each time step for Fi and T sm over 1000 s. 
The time average and maximum values of these errors are reported in Tab. H] 



Tab. 4: Mean and maximum values of the relative error of Fi and T sm for a duration of 
1000 s. . 



mesh 


e(Fi) 


&max 


c(2~' sm ) 


&max (Tam) 


mi 


0.11% 


0.75% 


0.68% 


4.34% 


m 2 


0.06% 


0.40% 


0.09% 


1.43% 



It is evident from Tabs. |3J and H] that the results given by mesh m 2 are sufficiently 
close to those given by mesh m 3 , which is two times finer. Thus, the size of mesh m 2 was 
adopted to perform the computations presented in this work, and all the meshes generated 
for the five enclosures have approximately 11,000 cells. 



5 Results 



The computations for the five enclosures were conducted assuming the same conditions 
and using identical numerical parameters. The time step size was fixed to a relatively 
small value At = 10~ 2 s, and at each time step, the convergence criteria of the SIMPLE 
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algorithm were fixed at 10~ 4 for the velocity residuals and 10~ 6 for the temperature 
residual. The latter criterion of convergence was fixed to a lower value in order to achieve 
a higher level of energy conservation in its both sensible and latent forms as they are of 
different orders of magnitude, furthermore, the released energy is tightly related to the 
amount of melted PCM. The residuals are calculated from Eq. f|T6|) as follows: 



where n c is the number of cells in the computational domain. Within each SIMPLE 
iteration, the liquid fraction update algorithm is stopped when the difference between 
two successive / values is less than 10~ 4 . At the first SIMPLE iteration of a time step, a 
few iterations are needed for / to converge (less than 3), and in the remaining SIMPLE 
iterations, only one iteration is needed. For numerical stability reasons, the imposed heat 
flux rate at the surface Sh is progressively raised during the first 5 seconds to reach its 
constant value q" = 10 4 W.m~ 2 . This practice enables easier convergence for the early 
time steps. 

5.1 Velocity field and melting front evolution 

We show the flow patterns represented by velocity vectors in the five studied enclosures in 
Fig. El This figure also gives the position of the melting front (or the solid phase zones) 
at five different instants during the process. At t = 50 s, a thin vertical layer of PCM has 
melted, and it has roughly the same size and the same form in all the enclosures. The 
flow field has a regular form with ascendant and descendant parallel fluid currents. Up 
to this time, the five enclosures have the same behavior, and the geometry does not play 
a role. At t — 100 s, more PCM has melted, especially near the upper part of the surface 
Sh- The behavior of the enclosures then differs visibly: the enclosures (a) and (b) that 
do not exceed the surface Sh have a melting front that extends more horizontally due to 
the presence of the upper boundary. Thus, the shapes of melted zones are comparable 
between enclosures (a) and (b) and among (c), (d) and (e). 

At t = 200 s, the melted PCM zone in the enclosure (e) extends upward with a 
round shape, and as a consequence, more solid PCM is exposed to the flow current of the 
liquid for this enclosure than the others (the liquid-solid interface is the most extended 
for enclosure (e)). From this instant on, the melting fronts in the five enclosures present 
an inclined curve because the melting is more advanced in the upper regions due to the 
thermal stratification (see Fig. |SJ). The flow patterns in the enclosure (e) show the 
formation of recirculation above the surface Sh- This recirculation, visible at t — 200 s 
(Fig. ED, is in fact an unsteady oscillatory phenomenon with a period less than 2 s, which 
is why it is not visible in the figure when t = 300 s. These oscillations are induced by 
the formation of a large liquid space above the surface Sh- At the upper extremity of 
Sh, the hot flow stream rising along the surface oscillates between two orientations: a 
vertical orientation, as observed at t = 300 s, and a roughly horizontal orientation, as 
observed at t = 200 s or t = 500 s, which results in a secondary recirculation in the 
upper-right corner. This unsteady behavior starts as early as t ~ 150 s for enclosure (e) 
but much later (t ~ 800 s) for enclosure (d) because the upper liquid space is narrower. 
The flows in the three other enclosures do not show any oscillatory characteristics. At 
t = 500 s, the melting in all the enclosures is almost complete. We can visually observe 




(21) 
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that enclosures (a) and (e) contain the least amount of solid PCM. This observation is 
confirmed quantitatively in Fig. [TTJ 

To give more quantitative details about the velocity field in the enclosures, we plot the 
vertical component of the velocity vector along the horizontal axis passing through the 
center of the surface Sh in Fig. [7] Initially (up to t — 100 s), the velocities are similar, 
so they are not shown. At t = 200 s, we observe substantial differences in the maximum 
value of the velocity near the boundary Sh, where enclosure (e) exhibits the highest 
value (6.5 x 10~ 3 m.s^ 1 ), and enclosure (a) exhibits the lowest value (4.7 x 10~ 3 m.s^ 1 ). 
The x positions of the positive maximum are nearly the same, whereas in the fluid zone 
near the melting front, the downward negative velocities have different profiles due to 
the differences in the position of the front. As the process continues (t = 300 s and 
t = 500 s), the velocity level decreases progressively, and a large zone with almost zero 
vertical velocity appears in the center of the enclosures. The maximum value in enclosure 
(e) at t = 500 s is 3.1 x 10 -3 m.s^ 1 . This velocity decrease is due to the reduction of 
the quantity of solid PCM, which is at the low temperature T = T m , in the enclosures 
compared to the initial amounts. As a result, a weaker temperature gradient exists in the 
liquid which results in weaker buoyancy forces. 

5.2 Temperature field in the enclosures 

Fig. [8] shows the temperature field at six different times during the melting process. As 
early as t = 200 s, we observe the establishment of thermal stratification in the melt, 
which is slightly perturbed in the upper region by the upward fluid currents and near 
the lateral boundaries by the convective heat exchange with the external media. Higher 
temperatures were observed in the liquid of enclosure (a), and lower temperatures were 
observed in enclosure (e). Quantitative changes in temperature are given in the next 
section. As for the melting front position, we observe similarities in the temperature fields 
between enclosures (a) and (b) and between enclosures (c) and (d), whereas enclosure 
(e) exhibits a distinct temperature pattern. 

5.3 Evolution of the parietal temperature 

The evolution of T smoy , the mean temperature of the surface Sh (Eq. ([20]) ). gives an 
overview of the impact of the flow and melting kinetics on the global temperature level 
(Fig. During the first minute, the mean temperature evolves identically and linearly 
for all the enclosures. This first phase is governed essentially by heat conduction; thus, 
it gives the same results for all five enclosures. After this phase of steep augmentation, 
Tsmoy continues to increase in enclosure (a) but with a weaker slope, but it decreases in all 
the other enclosures. The strongest decrease is observed for enclosure (e), and it lasts for 
approximately 250 s before increasing again. For enclosures (b), (c) and (d), this period 
is approximately 170 s. 

Another important temperature to monitor is the maximum temperature T max of the 
surface Sh because this temperature should be lower than the damage temperature of the 
device to be cooled. Fig. [10] plots the maximum temperature for the five enclosures. This 
maximum is always located at the top of Sh, as in Fig. [H] Its evolution shows differences 
between the enclosures as early as t = 30 s. For enclosure (a), T max exhibits two linear 
phases with positive slopes with inflexions at about t = 75 s and T max = 65°C. The 
evolution of T max for enclosure (b) has a horizontal plateau at 65°C, whereas for the 
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Fig. 6: Velocity vectors and melting fronts for the five enclosures after 50 and 100 s. 




Fig. 6: ( Cont.) Velocity vectors and melting fronts for the five enclosures after 200 and 300 s. 




Fig. 6: (Cont.) Velocity vectors and melting fronts for the five enclosures after 500 s. 
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Fig. 7: Vertical velocity components along the horizontal axis passing through the center of the 
surface Sh for the five enclosures after 200 s (a), 300 s (b) and 500 s (c). 
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other three enclosures, this plateau is 15°C lower and is slightly decreasing. The plateau 
for enclosure (e) (~ 200 s long) is the longest. 

The plateau in the T smoy and T max curves emerges because of the establishment of 
natural convection currents and the presence of the front of solid PCM. The solid PCM 
absorbs the thermal energy transported by the flow from the hot surface in the form of 
latent heat, and its melting releases cold liquid at T m , which increases the temperature 
gradient. The circular shape of enclosure (a) induces a larger distance between the solid 
PCM and the surface Sh- This distance is less important for enclosure (b), and its 
rectangular shape induces a plateau in the temperature profile and a lower value of T max . 
When the enclosures are flatter, as for (c) and (d), the results are better. The shape 
of the enclosure has another effect in terms of the extent of the melting front, which is 
more extended for the enclosures with the lowest aspect ratio. The more geometrically 
extended this melting front is in the enclosure, the lower and longer the T max plateau. By 
comparing the curves of enclosures (c) and (d) with (e), we can conclude that an efficient 
way to increase the extent of the front is to locate more PCM at the top of the enclosure 
where it is in contact with oscillatory natural convection streams. The higher temperature 
liquid is located in this zone because of thermal stratification, which increases the melting 
rate. 



5.4 Evolution of the global PCM melting 

The total liquid fraction F t defined by Eq. (TH?j) can be monitored to follow the evolution 
of the melting in the five enclosures, as shown in Fig. [TTJ We observe that the curve of 
Fi is linear and identical in all the enclosures until t « 150 s. After this time, the rate 
of melting decreases in all the enclosures except for enclosure (e), for which the curve 
inflexion occurs as late as t ~ 250 s. This enclosure has the fastest rate of melting, and 
the solid PCM completely disappears at t ~ 550 s. Enclosure (a) is also completely 
melted at this time, and even though it shows the highest maximum surface temperature 
(Fig. [10]), its mean surface temperature is lower than those of enclosures (b), (c) and (d) 
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after t ~ 350 s (Fig. [H]). Fig. 0can explain these differences. The bottom left corner of 
the rectangular enclosure (b) retains the PCM solid and limits its melting. In enclosures 

(c) and (d), more solid PCM is trapped in the lower part of the enclosure located below 
the surface Sh- This zone is less influenced by the hot liquid streams recirculating inside 
the enclosure, and thus, it is more difficult to melt. As a consequence, enclosures (c) and 

(d) are the last ones to completely melt, approximately 250 s after enclosures (a) and (e) 
are completely melted. 

5.5 Parietal heat transfer 

To quantify the heat transfer at the cooled surface, we define a global Nusselt number as 
follows: 

n" H 

Nu = q — (22) 

K\J- smoy ' rn ) 

Fig. [12] gives the evolution of the Nusselt number for the five enclosures during 
the melting process. After a fast transition period, where the Nu values are important 
{Tsmoy — T m « 0), and before they drop to a value of 20, we observe an increase toward 
a peak value followed by a progressive decrease toward an asymptotic value of 4. The 
characteristic "bumps" in the curves obtained for the different enclosures have different 
sizes and correspond to the coupled impact of the natural convection and latent heat of 
the melting of the PCM on the temperature of the cooled surface Sh- In agreement with 
the results of velocity (Fig. [7]) and parietal temperature (Fig. |9]), the most important 
Nu value is obtained for enclosure (e), and the lowest is obtained for enclosure (a) up to 
t ~ 350 s. After t ~ 350 s, enclosure (a) performs better than enclosures (b), (c) and 
(d) when the rate of fusion in these enclosures slows down. However, the most interesting 
behaviors expected from this mode of passive cooling are the fast reaction and limited 
maximum temperature. With regard to these two features, enclosure (a) is the worst 
choice and enclosure (e) is the best choice. We can also notice that for all the efficiency 
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Fig. 11: Evolution of the total liquid fraction in the five enclosures. 



indicators (T smoy , T max , Fi and Nu), the oblong enclosure (c) performs slightly better 
than the rectangular enclosure (d). 

6 Conclusions 

In this study, we numerically modeled the natural convection-dominated melting of a 
PCM inside an enclosure used to control the temperature of a surface releasing a heat 
flux. We investigated the impact of the shape of this enclosure and the relative position of 
the cooled surface on the flow and heat transfer. We developed a numerical model using 
a fixed grid enthalpy-porosity technique coupled with a three-dimensional flow solver 
based on an unstructured finite-volume method involving second-order accurate spatial 
and temporal numerical schemes. 

Five geometries containing the same quantity of PCM were compared. The geometries 
were rounded or rectangular, thick or thin, centered relative to the cooled surface or 
shifted upward vertically. The unsteady behaviors of the five enclosures were analyzed 
by examining the evolution of the liquid-solid interface and considering the forms and 
strengths of the natural convection flow. The fluid velocity increases for thinner enclosures 
that have a space above the cooled surface; thus, in the enclosure shifted upward (enclosure 
(e)) the maximum velocity is achieved. In spite of the presence of this flow, strong thermal 
stratification exists in the melt, which causes a hot spot located at the upper extremity 
of the surface to appear. 

Several efficiency indicators were monitored, including the mean and maximum tem- 
peratures, the total liquid fraction and the parietal Nusselt number. All of these indicators 
demonstrate the importance of the effect of the geometry on cooling efficiency. If we focus 
on the maximum temperature indicator, we see that at as early as t = 100 s, the differ- 
ence between the best and worst geometry choices ((e) and (a), respectively) is equal to 
16.5 °C, and this difference reaches 40 °C when t = 300 s, which is of great importance 
for the thermal protection of the heat-dissipating devices. We can summarize the most 
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important findings of this work as follows: 

• In thin enclosures, liquid PCM zones expand upward, while for wider enclosures, 
this expansion tends to be horizontal. Consequently, for the former enclosures, the 
solid cold PCM is in contact with ascending hot liquid streams, which contributes 
to the lowering of the global temperature. 

• The last zones to melt in the enclosures are located at the bottom and on the side 
opposite to the cooled surface. Thus, PCM enclosures should be designed without 
corners and without zones located below this surface. 

• A portion of the PCM should be placed above the cooled surface. 

• The use of rounded corners has a slight positive effect on the efficiency (e.g., Nu 
and T max ). 

Based on these findings, more efficient geometries can be designed; however, in the 
future, the impact of the geometry choice on the inverse process, i.e., solidification, should 
be taken into account and analyzed. Furthermore, the relations between the fraction 
of nanoparticles introduced and the thermophysical properties of the PCM such as the 
viscosity of the liquid phase, should be modeled and their impacts on hydrodynamics and 
heat transfer analyzed. These aspects constitute some of the goals of a future work. 
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A Test case of natural convection between a cylinder and a square 
duct 

Here we use the benchmark test case proposed by Demirdzic et al. pTl] to validate the 
implementation of natural convection in the code. A cylinder is placed inside a square 
enclosure of dimensions L x L with L = 1. The cylinder has a diameter of d/L = 0.4 and 
its center is displaced vertically by 5/L = 0.1. The cylinder is heated at Th = 1, while 
the vertical walls of the duct are cooled at T c = and the horizontal walls are assumed 



A Test case of natural convection between a cylinder and a square duct 



25 



to be adiabatic. The Rayleigh number based on L is 10 6 , and the Prandtl number is 
0.1. The benchmark solution of Demirdzic et al. [H] was obtained using a structured 
mesh of 256 x 128. Because of the symmetry of the problem, only one-half of the domain 
needs to be modeled. For our computation we used a mesh composed of 15, 366 cells, 
with quadrilateral cells around the walls and triangular cells elsewhere as showed in Fig. 
IT3l In this figure are presented the calculated streamlines and temperature contours that 
compare well with those of [41] . The local Nusselt number along the vertical wall of the 
duct is presented in Fig. [TH A very good agreement is found between our results and 
those of the benchmark. Furthermore, the calculated mean Nusselt number at this wall is 
6.7195 while the benchmark value is 6.7303 which gives a relative difference of less than 
0.1 %. 




Fig. 13: Calculated results of the natural convec- 
tion test case of Demirdzic et al. [4Tj . 
Right half: the used grid and tempera- 
ture contours (values from 0.05 to 0.95 
with a step of 0.1). Left half: streamline 
pattern. 



Fig. 14: Distribution of the local 
Nusselt number along the 
vertical cold wall compared 
to the benchmark data of 
Demirdzic et al. 



